
#### PACKAGES ####

library(tidyverse)
library(plm)
library(lme4)
library(clubSandwich)
library(modelsummary)

#### DATA GENERATION ####

conflictmerged <- readRDS("./conflictmerged_impute20.rds") %>%
  distinct(.keep_all = TRUE) %>%
  mutate(ics_number_main = as.factor(ics_number_main))

#### DESCRIPTIVE ####
desc1 <-  conflictmerged %>%
  select("Negotiation time" = diff, "Number of pages" = pages, "Productivity" = productivity_raw, "Productivity logged" = productivity, 
         "Flesch reading ease score" = re, "Edition" = edition, "ICS number" = ics_number_main,
         "Average distance between capitals" = distcap, "Share preferential trade agreements" = percent_pta, 
         "Avergae bilateral trade" = total_trade, "Average bilateral tariffs" = total_tariff,
         "Average regime distance" = ideological_diff, "Average UN voting distance" = un_diff, "Average UNGA mentions" = unga_speech, 
         "Share defensive alliances" = percent_alliance, "Share strategic rivals" = percent_rivalry) %>%
    as.data.frame() 
stargazer::stargazer(desc1, type = "text")

## ICC 

lmer(productivity ~ ideological_diff + (1 | committee), data = conflictmerged %>% 
       filter(year <= 2021) %>%
       mutate(ideological_diff = scale(ideological_diff, scale = TRUE, center = TRUE),
              re = scale(re, scale = TRUE, center = TRUE),
              edition = scale(edition, scale = TRUE, center = TRUE),
              productivity = scale(productivity, scale = TRUE, center = TRUE)))

0.4766/(0.4766 + 0.8078)  
# 0.3710682

# Committee explains about 37 percent of variation of productivity

##################################################### MULTILEVEL #####################################################

#### ECONOMIC ####

m1 <- lmer(productivity ~ distcap + factor(year) +
             edition + re + ics_number_main +
             (1 + distcap | committee), 
           data = conflictmerged %>% 
             filter(year <= 2020) %>%
             mutate(distcap = scale(distcap, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE) 
conflictmerged_scaled_1 <- conflictmerged %>% 
  select(productivity, distcap, year, committee,
         edition, ics_number_main, re) %>%
  filter(year <= 2020) %>%
  mutate(distcap = scale(distcap, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na() 

m2 <- lmer(productivity ~ percent_pta + factor(year) +
             edition + re + ics_number_main +
             (1 + percent_pta | committee), 
           data = conflictmerged %>% 
             filter(year <= 2020) %>%
             mutate(percent_pta = scale(percent_pta, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE)  
conflictmerged_scaled_2 <- conflictmerged %>% 
  select(productivity, percent_pta, year, committee,
         edition, ics_number_main, re) %>%
  filter(year <= 2020) %>%
  mutate(percent_pta = scale(percent_pta, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

m3 <- lmer(productivity ~ trade + as.factor(year) +
             edition + re + ics_number_main +
             (1 + trade | committee), 
           data = conflictmerged %>% 
             filter(year <= 2021) %>%
             mutate(trade = scale(total_trade, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE) 
conflictmerged_scaled_3 <- conflictmerged %>% 
  select(productivity, total_trade, committee, year,
         edition, ics_number_main, re) %>%
  filter(year <= 2021) %>%
  mutate(trade = scale(total_trade, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

m4 <- lmer(productivity ~ tariff + factor(year) +
             edition + re  + ics_number_main +
             (1 + tariff | committee), 
           data = conflictmerged %>% 
             filter(year <= 2021) %>%
             mutate(tariff = scale(total_tariff, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE) 
conflictmerged_scaled_4 <- conflictmerged %>% 
  select(productivity, total_tariff, committee, year,
         edition, ics_number_main, re) %>%
  filter(year <= 2021) %>%
  mutate(tariff = scale(total_tariff, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

### POLITICAL ###

m5 <- lmer(productivity ~ ideological_diff + as.factor(year) +
             edition + re + ics_number_main + 
             (1 + ideological_diff | committee),
           data = conflictmerged %>% 
             filter(year <= 2021) %>%
             mutate(ideological_diff = scale(ideological_diff, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE) 
conflictmerged_scaled_5 <- conflictmerged %>% 
  filter(year <= 2021) %>% 
  select(productivity, ideological_diff, committee, year,
         edition, ics_number_main, re) %>% 
  filter(year <= 2021) %>%
  mutate(ideological_diff = scale(ideological_diff, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

m6 <- lmer(productivity ~ un_diff + factor(year) +
             edition + re  + ics_number_main + 
             (1 + un_diff | committee), 
           data = conflictmerged %>% 
             filter(year <= 2014) %>%
             mutate(un_diff = scale(un_diff, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = TRUE)
conflictmerged_scaled_6 <- conflictmerged %>%
  select(productivity, un_diff, committee, year,
         edition, ics_number_main, re, n, year, committee) %>%
  filter(year <= 2014) %>%
  mutate(un_diff = scale(un_diff, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

m7 <- lmer(productivity ~ unga_speech + factor(year) +
             edition + re + ics_number_main +
             (1 + unga_speech | committee), 
           data = conflictmerged %>% 
             filter(year <= 2021) %>%
             mutate(unga_speech = scale(unga_speech, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE)
conflictmerged_scaled_7 <- conflictmerged %>% 
  select(productivity, unga_speech, committee, year,
         edition, ics_number_main, re) %>%
  filter(year <= 2021) %>%
  mutate(unga_speech = scale(unga_speech, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

m8 <- lmer(productivity ~ percent_rivalry  + factor(year) +
             edition + re  + ics_number_main +
             (1 + percent_rivalry | committee), 
           data = conflictmerged %>% 
             filter(year <= 2010) %>%
             mutate(percent_rivalry = scale(percent_rivalry, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE) 
conflictmerged_scaled_8 <- conflictmerged %>% 
  select(productivity, percent_rivalry, committee, year,
         edition, ics_number_main, re, n, year, committee) %>%
  filter(year <= 2010) %>%
  mutate(percent_rivalry = scale(percent_rivalry, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

m9 <- lmer(productivity ~ percent_alliance + factor(year) +
             edition + re  + ics_number_main +
             (1 + percent_alliance | committee), 
           data = conflictmerged %>% 
             filter(year <= 2018) %>%
             mutate(percent_alliance = scale(percent_alliance, scale = TRUE, center = TRUE),
                    re = scale(re, scale = TRUE, center = TRUE),
                    edition = scale(edition, scale = TRUE, center = TRUE),
                    productivity = scale(productivity, scale = TRUE, center = TRUE)),
           REML = FALSE) 
conflictmerged_scaled_9 <- conflictmerged %>% 
  select(productivity, percent_alliance, committee, year,
         edition, ics_number_main, re) %>%
  filter(year <= 2018) %>%
  mutate(percent_alliance = scale(percent_alliance, scale = TRUE, center = TRUE),
         re = scale(re, scale = TRUE, center = TRUE),
         edition = scale(edition, scale = TRUE, center = TRUE),
         productivity = scale(productivity, scale = TRUE, center = TRUE)) %>%
  drop_na()

varcov1 <- vcovCR(m1, cluster = conflictmerged_scaled_1$committee, type = "CR0")
varcov2 <- vcovCR(m2, cluster = conflictmerged_scaled_2$committee, type = "CR0")
varcov3 <- vcovCR(m3, cluster = conflictmerged_scaled_3$committee, type = "CR0")
varcov4 <- vcovCR(m4, cluster = conflictmerged_scaled_4$committee, type = "CR0")
varcov5 <- vcovCR(m5, cluster = conflictmerged_scaled_5$committee, type = "CR0")
varcov6 <- vcovCR(m6, cluster = conflictmerged_scaled_6$committee, type = "CR0")
varcov7 <- vcovCR(m7, cluster = conflictmerged_scaled_7$committee, type = "CR0")
varcov8 <- vcovCR(m8, cluster = conflictmerged_scaled_8$committee, type = "CR0")
varcov9 <- vcovCR(m9, cluster = conflictmerged_scaled_9$committee, type = "CR0")

varcov1 <- as.matrix(varcov1)
varcov2 <- as.matrix(varcov2)
varcov3 <- as.matrix(varcov3)
varcov4 <- as.matrix(varcov4)
varcov5 <- as.matrix(varcov5)
varcov6 <- as.matrix(varcov6)
varcov7 <- as.matrix(varcov7)
varcov8 <- as.matrix(varcov8)
varcov9 <- as.matrix(varcov9)

cm <- c("distcap" = "Standardized coefficient",
        "percent_pta" = "Standardized coefficient",
        "count_pta" = "Standardized coefficient",
        "trade" = "Standardized coefficient",
        "tariff" = "Standardized coefficient",
        "ideological_diff" = "Standardized coefficient",
        "un_diff" = "Standardized coefficient",
        "unga_speech" = "Standardized coefficient",
        "percent_rivalry" = "Standardized coefficient",
        "percent_alliance" = "Standardized coefficient",
        "count_rivalry" = "Standardized coefficient",
        "count_alliance" = "Standardized coefficient",
        "count_co_patents" = "Standardized coefficient",
        "count_dom_own" = "Standardized coefficient",
        "count_for_own" = "Standardized coefficient",
        "co_patents" = "Standardized coefficient",
        "dom_own" = "Standardized coefficient",
        "n" = "Number of committee members",
        "standards" = "Standard production per year",
        "re" = "Complexity",
        "mean_chars" = "Abstract length",
        "pages" = "Pages",
        "statusWithdrawn" = "Status - Withdrawn",
        "edition" = "Edition")

models <- list(
  "Average geographic distance" = m1, 
  "Share PTAs" = m2, 
  "Average trade" = m3,
  "Average tariffs" = m4,
  "Average regime distance" = m5,
  "Average UN voting distance" = m6,
  "Average shared UNGA mentions" = m7,
  "Share strategic rivals" = m8,
  "Share defensive alliances" = m9
)

rows <- tibble::tribble(~start, ~m1, ~m2, ~m3, ~m4, ~m5, ~m6, ~m7, ~m8, ~m9,
                        "Time series", "2004-2020", "2004-2020",  "2004-2021", "2004-2021", "2004-2021", "2004-2014", "2004-2021", "2004-2010", "2004-2018") #"2004-2020", "2004-2020")

modelsummary(
  models,
  coef_map = cm,
  vcov = list(varcov1, varcov2, varcov3, varcov4, 
              varcov5, varcov6, varcov7, varcov8, varcov9),
  stars = TRUE,
  gof_omit = "AIC|BIC|Std.Errors",
  add_rows = rows,
  output = "default",
  notes = c("Fixed effects by year, random effects by committee.",
            "Clustered standard errors by committee.",
            "Includes standards that are Published and Withdrawn."))

# Interpreting coefficients (from standardized coef. and logged dep.var.)
sd_productivity <- sd(conflictmerged$productivity, na.rm = TRUE)
standardized_coefficient <- -0.027
unstandardized_coefficient <- standardized_coefficient * sd_productivity
percentage_change <- (exp(unstandardized_coefficient) - 1) * 100
print(c(unstandardized_coefficient, percentage_change))


fig1 <- broom.mixed::tidy(m1) %>%
  filter(term == "distcap") %>%
  bind_rows(broom.mixed::tidy(m2) %>%
              filter(term == "percent_pta")) %>%
  bind_rows(broom.mixed::tidy(m3) %>%
              filter(term == "trade")) %>%
  bind_rows(broom.mixed::tidy(m4) %>%
              filter(term == "tariff")) %>%
  bind_rows(broom.mixed::tidy(m5) %>%
              filter(term == "ideological_diff")) %>%
  bind_rows(broom.mixed::tidy(m6) %>%
              filter(term == "un_diff")) %>%
  bind_rows(broom.mixed::tidy(m7) %>%
              filter(term == "unga_speech")) %>%
  bind_rows(broom.mixed::tidy(m8) %>%
              filter(term == "percent_rivalry")) %>%
  bind_rows(broom.mixed::tidy(m9) %>%
              filter(term == "percent_alliance")) %>%
  mutate(conf.low = c(conf_int(m1, vcov = "CR0", cluster = conflictmerged_scaled_1$committee, coefs = "distcap")[[5]],
                      conf_int(m2, vcov = "CR0", cluster = conflictmerged_scaled_2$committee, coefs = "percent_pta")[[5]],
                      conf_int(m3, vcov = "CR0", cluster = conflictmerged_scaled_3$committee, coefs = "trade")[[5]],
                      conf_int(m4, vcov = "CR0", cluster = conflictmerged_scaled_4$committee, coefs = "tariff")[[5]],
                      conf_int(m5, vcov = "CR0", cluster = conflictmerged_scaled_5$committee, coefs = "ideological_diff")[[5]],
                      conf_int(m6, vcov = "CR0", cluster = conflictmerged_scaled_6$committee, coefs = "un_diff")[[5]],
                      conf_int(m7, vcov = "CR0", cluster = conflictmerged_scaled_7$committee, coefs = "unga_speech")[[5]],
                      conf_int(m8, vcov = "CR0", cluster = conflictmerged_scaled_8$committee, coefs = "percent_rivalry")[[5]],
                      conf_int(m9, vcov = "CR0", cluster = conflictmerged_scaled_9$committee, coefs = "percent_alliance")[[5]])) %>%
  mutate(conf.high = c(conf_int(m1, vcov = "CR0", cluster = conflictmerged_scaled_1$committee, coefs = "distcap")[[6]],
                       conf_int(m2, vcov = "CR0", cluster = conflictmerged_scaled_2$committee, coefs = "percent_pta")[[6]],
                       conf_int(m3, vcov = "CR0", cluster = conflictmerged_scaled_3$committee, coefs = "trade")[[6]],
                       conf_int(m4, vcov = "CR0", cluster = conflictmerged_scaled_4$committee, coefs = "tariff")[[6]],
                       conf_int(m5, vcov = "CR0", cluster = conflictmerged_scaled_5$committee, coefs = "ideological_diff")[[6]],
                       conf_int(m6, vcov = "CR0", cluster = conflictmerged_scaled_6$committee, coefs = "un_diff")[[6]],
                       conf_int(m7, vcov = "CR0", cluster = conflictmerged_scaled_7$committee, coefs = "unga_speech")[[6]],
                       conf_int(m8, vcov = "CR0", cluster = conflictmerged_scaled_8$committee, coefs = "percent_rivalry")[[6]],
                       conf_int(m9, vcov = "CR0", cluster = conflictmerged_scaled_9$committee, coefs = "percent_alliance")[[6]]))

figure1 <- fig1 %>%
  mutate(group = ifelse(term %in% c("ideological_diff", "un_diff", "unga_speech", "percent_rivalry", "percent_alliance"), "Political",
                        ifelse(term %in% c("trade", "tariff", "distcap", "percent_pta"), "Economic", 
                               ifelse(term %in% c("co_patents", "dom_own"), "Scientific",
                                      ifelse(term %in% c("pdi_diff", "mas_diff", "uai_diff", "idv_diff", "ltowvs_diff", "ivr_diff"), "Cultural",
                                      "none"))))) %>%
  mutate(term = ifelse(term == "ideological_diff", "Average regime distance",
                       ifelse(term == "un_diff", "Average UN voting distance",
                              ifelse(term == "unga_speech", "Average UNGA mentions",
                                     ifelse(term == "percent_rivalry", "Share strategic rivals",
                                            ifelse(term == "percent_alliance", "Share defensive alliances",
                                                   ifelse(term == "distcap", "Average geographic distance",
                                                          ifelse(term == "percent_pta", "Share preferential trade agreements",
                                                                 ifelse(term == "trade", "Average bilateral trade",
                                                                        ifelse(term == "tariff", "Average bilateral tariffs",
                                                                        ifelse(term == "co_patents", "Patent collaborations",
                                                                               ifelse(term == "dom_own", "Co-owned inventions",
                                                                                      ifelse(term == "pdi_diff", "Power Distance Index",
                                                                                             ifelse(term == "mas_diff", "Masculinity vs. Femininity",
                                                                                                    ifelse(term == "uai_diff", "Uncertainty Avoidance Index",
                                                                                                           ifelse(term == "idv_diff", "Individualism vs. Collectivism",
                                                                                                                  ifelse(term == "ltowvs_diff", "Long-Term vs. Short-Term Orientation",
                                                                                                                         ifelse(term == "ivr_diff", "Indulgence vs. Restraint", 
                                                                                                                                NA)))))))))))))))))) %>%
  ggplot(aes(term, estimate, label = group)) + 
  geom_point(size = 4) + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), colour = "gray", linewidth = 1.5) + 
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 1.2) +
  labs(x = "", y = "Committee producitvity")  + 
  facet_wrap(~group, scales ="free_y", nrow = 3) +
  coord_flip() +
  theme_bw() +
  theme(axis.title=element_text(size=20, face="bold"),
        axis.text =element_text(size=rel(3)),
        strip.text.x = element_text(size = 30)) 

figure1


